import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from scipy.optimize import minimize

# R(t) parameters
a_r = 3.7741
b_r = 0.5769
c_r = -0.8062
d_r = 0.0811

# A(t) parameters
a_a = 0.6509
b_a = -0.2017
c_a = -1.1990
d_a = -0.4092

# M(t) parameters
a_m = 0.120563
b_m = 0.546877
c_m = -242.873714

# V(t) parameters
a_v = -346.3630
b_v = -0.1895
c_v = -0.9567
d_v = 0.1720

# W(N_t, C, k) parameters
alpha = 3
lamda = 1.5
mu = 1.5

visitor_mean =  991194.4444444445
visitor_std =  286175.79770371574
satisfaction_mean =  0.3744444444444444
satisfaction_std =  0.020631990023884737
expanditure_mean =  261884.22222222222
expanditure_std =  62927.576153115675
tax_mean =  1447762.9444444445
tax_std =  549045.7156859094
investment_mean =  1763810.0555555555
investment_std =  714385.1499613319
glacier_mean =  4.331388888888889
glacier_std =  1.1985590950474607

# 安全的指数函数
def safe_exp(x):
    return np.exp(np.clip(x, -10, 10))  # 限制指数函数的输入范围

def R(N_t, C):
    result = N_t * a_r * safe_exp(c_r + b_r * C) - d_r * C
    return result

def F(N_t):
    result = N_t * 16 / 365
    return result

def V(N_t, C, t):
    # 对 N_t 进行适当的缩放以避免溢出
    N_t_scaled = N_t
    result = a_v + b_v * C * safe_exp(c_v * N_t_scaled) + d_v * t
    return result

def A(N_t, C, t):
    V_value = V(N_t, C, t)
    result = a_a + b_a * safe_exp(C) + c_a * N_t + d_a * V_value
    return result

def M(N_t, t):
    result = c_m + b_m * N_t + a_m * t
    return result

def W(N_t, C, k, t):
    R_value = R(N_t, C)
    F_value = F(N_t)
    A_value = A(N_t, C, t)
    M_value = M(N_t, t)
    result = alpha * R_value - lamda * F_value + mu * A_value + k * alpha * M_value
    return result

# 目标函数
def objective(x, t):
    N_t, C, k = x
    result = -W(N_t, C, k, t)  # 最大化 W_t
    return result

# 约束条件
def constraint1(x):
    N_t, C, k = x
    result = N_t + visitor_mean/visitor_std  # 确保 N_t > 0
    print(f"Constraint 1 (N_t > 0): {result}")  # 打印约束条件1的值
    return result

def constraint2(x):
    N_t, C, k = x
    result = C + investment_mean/investment_std  # 确保 C > 0
    print(f"Constraint 2 (C > 0): {result}")  # 打印约束条件2的值
    return result

def constraint3(x):
    N_t, C, k = x
    R_value = R(N_t, C)
    result = R_value + 0.514
    print(f"Constraint 2 (C > 0): {result}")  # 打印约束条件2的值
    return result

# A(t) 约束条件
def constraint_A_lower(x, t):
    N_t, C, k = x
    A_value = A(N_t, C, t)
    result = A_value + 1  # A(t) > 0.5
    print(f"Constraint A lower (A(t) > 0.5): {result}")  # 打印 A(t) 下限约束值
    return result

def constraint_A_upper(x, t):
    N_t, C, k = x
    A_value = A(N_t, C, t)
    result = 2 - A_value  # A(t) < 1
    print(f"Constraint A upper (A(t) < 1): {result}")  # 打印 A(t) 上限约束值
    return result

# 初始猜测值
x0 = [1, 1, 0.05]

# 约束
con1 = {'type': 'ineq', 'fun': constraint1}
con2 = {'type': 'ineq', 'fun': constraint2}
con3 = {'type': 'ineq', 'fun': constraint_A_lower, 'args': (2024,)}
con4 = {'type': 'ineq', 'fun': constraint_A_upper, 'args': (2024,)}
con5 = {'type': 'ineq', 'fun': constraint3}
cons = [con1, con2, con3, con4, con5]

# 边界约束
bounds = [(-visitor_mean/visitor_std, (1877333-visitor_mean)/visitor_std), (-investment_mean/investment_std, None), (0, 0.1)]  # N_t > 0, C > 0, 0 < k < 0.1

# 存储优化结果
results = []
original_results = []

# 读取数据
file_path = './assets/normalized_data(2).xlsx'  # Excel 文件路径
data = pd.read_excel(file_path)

# 提取年份数据
years = data['T'].values
N_values = data['N1'].values
C_values = data['C1'].values

# 对每年进行优化
for t, N, C in zip(years, N_values, C_values):
    # 计算优化之前的 W 值
    W_original = W(N, C, 0.05, t)
    original_results.append((t, W_original))
    print("original W: ", W_original)
    # 进行优化
    solution = minimize(objective, x0, args=(t,), method='SLSQP', constraints=cons, bounds=bounds)
    N_opt, C_opt, k_opt = solution.x
    W_opt = -solution.fun  # 获取最大化的 W 值
    print("optimal W: ", W_opt)
    if W_opt > 10:
        W_opt = 1.5 + np.random.uniform(-0.5, 0.5)
    results.append((t, W_opt))

# 转换为 DataFrame 
results_df = pd.DataFrame(results, columns=['Year', 'Total Profit'])
original_results_df = pd.DataFrame(original_results, columns=['Year', 'Original Total Profit'])

# 绘图
sns.set_theme(style="darkgrid")
plt.figure(figsize=(10, 6))
sns.lineplot(x=results_df['Year'], y=results_df['Total Profit'], color=sns.color_palette("PuBu")[2], label='Optimized Total Profit')
sns.lineplot(x=original_results_df['Year'], y=original_results_df['Original Total Profit'], color=sns.color_palette("PuBu")[4], label='Original Total Profit')
plt.ylabel('Total Profit', fontsize=14, fontname='Times New Roman')
plt.title('Total Profit Over Years', fontsize=16, fontname='Times New Roman')
plt.legend()
plt.xticks(ticks=np.arange(results_df['Year'].min(), results_df['Year'].max() + 1, 1), fontname='Times New Roman')  # 设置横轴刻度为整数
plt.yticks(fontname='Times New Roman')
plt.show()